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Abstract 

We propose new methods for the numerical continuation of point-to-cycle connecting orbits 
in 3-dimensional autonomous ODE's using projection boundary conditions. In our approach, 
the projection boundary conditions near the cycle are formulated using an eigenfunction of the 
associated adjoint variational equation, avoiding costly and numerically unstable computations 
of the monodromy matrix. The equations for the eigenfunction are included in the defining 
boundary- value problem, allowing a straightforward implementation in AUTO, in which only 
the standard features of the software are employed. Homotopy methods to find connecting 
orbits are discussed in general and illustrated with several examples, including the Lorenz equa- 
tions. Complete auto demos, which can be easily adapted to any autonomous 3-dimensional 
ODE system, are freely available. 

Keywords: boundary value problems, projection boundary conditions, point-to-cycle connec- 
tions, global bifurcations 
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1 Introduction 

Many interesting phenomena in ODE systems can 
only be understood by analyzing global bifurca- 
tions. Examples of such are the occurrence and 
disappearance of chaotic behaviour. For exam- 
ple, the classical Lorenz attractor appears in a 
sequence of bifurcations, where homoclinic orbits 
connecting a saddle equilibrium to itself and het- 
eroclinic orbits connecting an equilibrium point 
with a saddle cycle, are involved (Afraimovich et 
al., 1977). In the ecological context, Boer et al. 
(1999, 2001) showed that regions of chaotic be- 
haviour in parameter space in some food chain 
models are bounded by bifurcations of point-to- 
cycle and cycle-to-cycle connections. 

Thus, in order to gain more knowledge about 
the global bifurcation structure of a model, in- 
formation is required on the existence of homo- 
clinic and heteroclinic connections between equi- 
libria and/or periodic cycles. The first type is a 
connection that links an equilibrium or a cycle to 
itself (asymptotically bi-stable, so it necessarily 
has nontrivial stable and unstable invariant man- 
ifolds). The second type is a connection that links 
an equilibrium or a cycle to another equilibrium 
or cycle. 

The continuation of connecting orbits in ODE 
systems has been notoriously difficult. Doedel 
and Friedman (1989) and Beyn (1990) developed 
direct numerical methods for the computation of 
orbits connecting equilibrium points and their as- 
sociated parameter values, based on truncated 
boundary value problems with projection bound- 
ary conditions. Moreover, Doedel, Friedman and 
Monteiro (1993) have proposed efficient methods 
to find starting solutions by successive continua- 
tions (homotopies) . These continuation methods 
have been implemented in HomCont, as incor- 
porated in AUTO (Doedel et al., 1997; Champneys 
and Kuznetsov, 1994; Champneys et al., 1996). 
HomCont is only suitable for the continuation of 
homoclinic point-to-point and heteroclinic point- 
to-point connections. 

More recently, significant progress has been 



made in the continuation of homoclinic and het- 
eroclinic connections involving cycles. Dieci and 
Rebaza (2004) developed a method based on ear- 
lier works by Beyn (1994) and Pampel (2001). 
Their method is also based on projection bound- 
ary conditions, but uses an ad hoc multiple shoot- 
ing technique and requires the numerical determi- 
nation of the monodromy matrix associated with 
the periodic cycles involved in the connection. 

In this paper, we propose new methods for the 
numerical continuation of point-to-cycle connec- 
tions in 3-dimensional autonomous ODE's using 
projection boundary conditions. In our approach, 
the projection boundary conditions near each cy- 
cle are formulated using an eigenfunction of the 
associated adjoint variational equation, avoiding 
costly and numerically unstable computation of 
the monodromy matrix. Instead, the equations 
for the eigenfunction are included in the defining 
boundary-value problem, allowing a straightfor- 
ward implementation in AUTO. 

This paper is organized as follows. In Sec- 
tion [2] we recall basic properties of the projection 
boundary condition method to continue point-to- 
cycle connections. In Section [3] this method is 
adapted to efficient numerical implementation in 
a special - but important - 3D case. Homotopy 
methods to find connecting orbits are discussed 
in Section HI Section [5] demonstrates that the 
algorithms allow for a straightforward implemen- 
tation in AUTO, using only the basic features of 
this software. Three well-known examples (the 
three-dimensional Lorenz system, the electronic 
circuit model of Freire et al., 1993, and the stan- 
dard three-level food chain model based on the 
Rosenzweig-MacArthur (1963) system) are used 
in Section [H] to illustrate the power of the new 
methods. 

This is Part I of a sequel of two papers. Part 
II will deal with cycle-to-cycle connections in 3D 
systems. 
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Figure 1: Point-to-cycle connecting orbits in M : (a) ra^ = 1; (b) n„ = 2 



2 Truncated BVP's with pro- 
jection BC's 

Before presenting a BVP for a point-to-cycle con- 
nection, we set up some notation. 

Consider a general system of ODE's 



du 
Itt 



fiu,a), 



(1) 



where / : x ^ is a sufficiently smooth 
function of the state variables u G M"' and the 
control parameters a G M^. Denote by cp^ the 
(local) flow generated by ([T])0. 

Let 0~ be either a saddle or a saddle- focus 
equilibrium, say ^, and let be a hyperbolic 
saddle limit cycle of ([1]). A solution u(t) of ([1]) 
defines a connecting orbit from 0~ to if 



lim dist(M(t),0^) = 

t— »±oo 



(2) 



(see Figure [T] for illustrations). Since u{t + r) 
satisfies ([1]) and ([2]) for any phase shift r, an ad- 
ditional scalar phase condition 



ip[u, a] = 



(3) 



^Whenever possible, we will not indicate explicitly the 
dependence of various objects on system parameters. 



is needed to ensure uniqueness of the connecting 
orbit. This condition will be specified later. 

For numerical approximation, the asymptotic 
conditions ([2]) are replaced by projection boundary 
conditions at the end-points of a large truncation 
interval [r„,r_|_]: The points u{t_) and u{t^) are 
required to belong to the linear subspaces that 
are tangent to the unstable and stable invariant 
manifolds of 0~ and 0~^, respectively. 

Let n~ be the dimension of the unstable in- 
variant manifold W!^ of ^, i.e., the number of 
eigenvalues A~ of the Jacobian matrix = D^f 
evaluated at the equilibrium which satisfy 

3?(A-) > 0. 

Denote by x^{t) a periodic solution (with min- 
imal period T+) corresponding to and intro- 
duce the monodromy matrix 



T+, 



i.e., the linearization matrix of the T"''-shift along 
orbits of ([1]) at point Xq = x~^{0) G O"*". Its eigen- 
values /i"*" are called the Floquet multipliers; ex- 
actly one of them equals 1, due to the assumption 
of hyperbolicity. Let = -|- 1 be the dimen- 
sion of the stable invariant manifold Wl of the 
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cycle O"*"; here is the number of its multiphers 
satisfying 

l/i+l < 1. 

A necessary condition to have an isolated family 
of point-to-cycle connecting orbits of ([T]) is that 
(see Beyn (1994)) 

p = n — mj" — n~ + 2 (4) 

The projection boundary conditions in this 
case can be written as 

L-(n(r_)-O=0, (5a) 
L+(M(r+) -x+(0)) = , (5b) 

where L~ is a {n—n~)xn matrix whose rows form 
a basis in the orthogonal complement of the linear 
subspace that is tangent to at ^. Similarly, 
is a (n — m+) x n matrix, such that its rows 
form a basis in the orthogonal complement to the 
linear subspace that is tangent to of O"^ at 
x+(0). 

It can be proved that, generically, the trun- 
cated BVP composed of ([T]), a truncation of ([3]), 
and has a unique solution family {it, a), pro- 
vided that has a connecting solution family 
satisfying ([3]) and (jlj). 

The truncation to the finite interval [r_,r4.] 
implies an error. If m is a generic connecting so- 
lution to ([T]) at parameter value a, then the fol- 
lowing estimate holds: 

||(n|[._,.^],a) - {u,a)\\ < Ce"2-in(;.-|r_|,;.+ |.^|)^ 

where || ■ || is an appropriate norm in the space 
C^([r_, r+], M") X W, M|[r_,r+] is the restriction of 
u to the truncation interval, and fi± are deter- 
mined by the eigenvalues of the Jacobian matrix 
and the monodromy matrix. See Pampel (2001) 
and Dieci and Rebaza (2004) for exact formula- 
tions, proofs, and references to earlier contribu- 
tions. 

3 New defining systems in 

Here we explain how the projection boundary con- 
ditions ([5]) can be implemented efficiently in a 



special - but important - case n = 3. Thereafter 
we specify the defining system used to continue 
connecting orbits in 3D-0DE example systems 
with AUTO. A saddle cycle 0+ in such systems 
always has m+ = = 2. 

3.1 The equilibrium-related part 

The equilibrium point ^, an appropriate solution 
of /(^, a) = 0, cannot be found by time-integration 
methods because it is a saddle. There are two dif- 
ferent types of saddle equilibria that can be con- 
nected to saddle cycles in 3D-0DE's. These are 
distinguished by the dimension n~ of the unsta- 
ble invariant manifold of ^: We have either 
n~ = 1 OT n~ = 2 (see Figured]). In the former 
case, the connection is structurally unstable (has 
codim 1) and, according to (jlj), we need two free 
system parameters for its continuation {p = 2). In 
the latter case, however, the connection is struc- 
turally stable and can be continued, generically, 
with one system parameter [p = 1). There is 
a small difference in the implementation of the 
projection boundary condition fl5al) in these two 
cases. 

If ?T,~ = 1 (see Figure [2](a)), then the follow- 
ing explicit projection boundary condition replaces 

u{T.)=^ + ev, (6) 

where e > is a given small number, and v G M.^ 
is a unit vector that is tangent to W!^ at C,- Notice 
that this fixes the phase of the connecting solution 
u, so that ([3]) becomes (l5a|) in this case. The vec- 
tor f in ([6]) is, of course, a normalized eigenvector 
associated with the unstable eigenvalue A„ > of 
the Jacobian matrix evaluated at the equilib- 
rium. Hence, we can use the following algebraic 
system to continue ^, v and A„ simultaneously: 

/(e,«) = 0, 

f^{^,a)v-X^v = 0, (7) 
{v,v)-l = 0, 

where {x, u) = x^u is the standard scalar product 
in W. 
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(a) 



(b) 



Figure 2: EVP's to approximate connecting orbits: (a) = 1; (b) = 2. 



If n" = 2 (see Figure [2](b)), then H^" is or- 
thogonal to an eigenvector v of the transposed Ja- 
cobian matrix corresponding to its eigenvalue 
As < 0, so that (l5al) can be written as 

{v,u{t^)-O = 0. (8) 

To continue ^,v, and A^, we use a system similar 
to ([7]), namely: 

= 0, 

f^{^,a)v-X,v = 0, (9) 
{v,v)-l = 0. 

As a variant of the phase condition ([3]) in this 
case, we can use the linear condition 

(r/,n(r_)-O=0, (10) 

which places the starting point of the truncated 
connecting solution in a plane containing the equi- 
librium ^ and orthogonal to a fixed vector t] (not 
coUinear with v). 

3.2 The cycle and eigenfunctions 

The heteroclinic connection is linked on the other 
side to a saddle limit cycle O"*" (see Figure [2]). 



Thus, we also need a EVP to compute it. We use 
the standard periodic EVP: 

f x+ - /(x+,a) = , , . 

\ x+(0)-x+(T+) =0, ^ ' 

which is augmented by an appropriate phase con- 
dition that makes its solution unique. This phase 
condition is actually a boundary condition for the 
truncated connecting solution, and will be intro- 
duced below. 

To set up the projection boundary condition 
for the truncated connecting solution u near O"*", 
we also need a vector, say w{0), that is orthogonal 
at x(0) to the stable manifold of the saddle 
limit cycle 0+ (see Figure [2]). It is well known 
that w{0) can be obtained from an eigenf unction 
w{t) of the adjoint variational problem associated 
with (ITTl) . corresponding to its eigenvalue 

1 

where /i^ is a multiplier of the monodromy matrix 
satisfying 

(see Appendix). The corresponding EVP is 

w + fj{x+,a)w = 0, 
w(T+) - /iw(0) = 0, (12) 
{wiO),wm-l = 0, 
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where is the solution of (ITTI) . In our implemen- 
tation the above BVP is replaced by an equivalent 
BVP 



or 



with 



w + {x~^ , a)w + Xw = 0, 

w{T+) - sw{0) = , 

{w{0),w{0))-l = 0, 

where s = sign /i = ±1 and 



A = In \fi\ 

(see Appendix). In (fT3l) . the boundary conditions 
become periodic or anti-periodic, depending on 
the sign of the multiplier fi, while the logarithm of 
its absolute value appears in the variational equa- 
tion. This ensures high numerical robustness. 

Given w satisfying f|T3l) . the projection bound- 
ary condition fl5bl) becomes 



(14) 



(ti;(0),n(r+)-x+(0)) = 0. 



3.3 The connection 

Finally, we need a phase condition to select a 
unique periodic solution among those which sat- 
isfy ( ITTI) . i.e., to fix a base point Xq = x~^{0) 
on the cycle O"*" (see Figure E]). Usually, an in- 
tegral condition is used to fix the phase of the 
periodic solution. For the point-to-cycle connec- 
tion, however, we need a new condition, since the 
end point near the cycle should vary freely. To 
this end we require the end point of the connec- 
tion to belong to a plane orthogonal to the vector 
Jq = f{x~^{0),a). This gives the following BVP 
for the connecting solution: 



u 



f{u,a) 
x+(0)) 



(/(x+(0),a),M(r+) 
3.4 The complete BVP 








(15) 



The complete truncated BVP to be solved numer- 
ically consists of ([7j), with 



uiO) = ^ + ev, 



(16) 



{v,u{0) 
(r/,M(0) 



= 
= 0, 



(17a) 
(17b) 



(13) as well as 





x+ -T+f{x+,a) 


= 0, 


(18a) 




x+(0) - x+(l) 


= 0, 


(18b) 


{w 


(0),m(1)-x+(0)) 


= 0, 


(18c) 




'f^{x^,a)w + Xw 


= 0, 


(18d) 




w{l) - sw{0) 


= 0, 


(18e) 




{w{0),w{0))-l 


= 0, 


(18f) 




ii — Tf{u, a) 


= 0, 


(18g) 




,«),m(1)-x+(0)) 


= 0. 


(18h) 



Here the time variable is scaled to the unit inter- 
val [0,1], so that both the cycle period and 
the connecting time T become parameters. 

If the connection time T is fixed at a large 
value, this BVP allows to continue simultaneously 
the equilibrium ^, its eigenvalue A^^ or A^, the cor- 
responding eigenvector v, the periodic solution x~^ 
corresponding to the limit cycle 0~^, its period 
T"*", the logarithm of the absolute value of the 
unstable multiplier of this cycle, the correspond- 
ing scaled eigenfunction w, as well as (a trunca- 
tion of) the connecting orbit u. These objects 
become functions of one system parameter (when 
dimW^ = 2) or two system parameters (when 
dimiy^ = 1). These free system parameters are 
denoted as Oj. 

If dim = 2 then, generically, limit points 
(folds) are encountered along the solution fam- 
ily. These can be detected, located accurately, 
and subsequently continued in two system pa- 
rameters, say, {ai,a2), using the standard fold- 
following facilities of auto. 

4 Starting strategies 

The BVP's specified above can only be used if 
good starting data are available. This can be 
problematic, since global objects - a saddle cycle 
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and a connecting orbit - are involved. However, a 
series of successive continuations in auto can be 
used to generate all necessary starting data, given 
little a priori knowledge about the existence and 
location of a heteroclinic point-to-cycle connec- 
tion. 

4.1 The equilibrium and the cycle 

The equilibrium ^, its unstable or stable eigen- 
value, as well as the corresponding eigenvector or 
adjoint eigenvector can be calculated using maple 
or MATLAB. Alternatively, this saddle equilib- 
rium can often be obtained via continuation of 
a stable equilibrium family through a limit point 
(fold) bifurcation. 

To obtain the limit cycle 0~^, one can continue 
numerically (with auto or content, for exam- 
ple) a limit cycle born at a Hopf bifurcation to an 
appropriate value of a, from where we start the 
successive continuation. 

4.2 Eigenfunctions 

In the first of such continuations, the periodic so- 
lution corresponding the limit cycle at the par- 
ticular parameter values is used to get an eigen- 
function. To explain the idea, let us begin with 
the original adjoint eigenfunction w. Consider the 
periodic EVP fll8ap - fll8bl) for the cycle, to which 
the standard integral phase condition is added, 

[\iUr),xHT)) = 0, (19) 
Jo 

as well clS db EVP similar to (fT2|) . namely: 

r w + T+f:^{x+,a)w = , 

I w{l)-fiw{0) = 0, (20) 

( {w{0),w{0)) -h = 0. 

In (fT9|) . x^i^ is a reference periodic solution, typ- 
ically the one in the preceding continuation step. 
The parameter h in (120!) is a homotopy parame- 
ter, that is set to zero initially. Then (1201) has a 
trivial solution 

w{t) = 0, h = 0, 



for any real This family of trivial solutions 
parametrized by fi can be continued in auto us- 
ing a EVP consisting of ffTTl) (with scaled time 
variable t), (HM . and fl2Ul) with free parameters 
(/i, h) and fixed a. A Floquet multiplier of the ad- 
joint system then corresponds to a branch point 
at /ii along this trivial solution family (see Ap- 
pendix). AUTO can accurately locate such a point 
and switch to the nontrivial branch that emanates 
from it. Continuing this secondary family in (/i, h) 
until, say, the value h = 1 is reached, gives a non- 
trivial eigenfunction w corresponding to the mul- 
tiplier fii. Note that in this continuation the value 
of fi remains constant, yU = /ii, up to numerical 
accuracy. 

The same method is applicable to obtain a 
nontrivial scaled adjoint eigenfunction. For this, 
the EVP 

( w + T+f^{x+,a)w + Xw = 0, 

<j w{l)-sw{0) = 0, (21) 

[ {wiO),w{0))-h = 0, 

where s = sign(yU,), replaces (!20|) . A branch point 
at Ai then corresponds to the adjoint multiplier 
se^^. Eranch switching then gives the desired 
eigendata. 

4.3 The connection 

Sometimes, an approximation of the connecting 
orbit can be obtained by time-integration of ([T]) 
with a starting point satisfying ([H]) or ([HD and flTU]) . 
These data (the periodic solution corresponding 
to the limit cycle, its nontrivial eigenfunction, 
and the integrated connecting orbit) must then 
be merged, using the same scaled time variable 
and mesh points. This only works for non-stiff 
systems provided that the connecting orbit and 
its corresponding parameter values are known a 
priori with high accuracy, which is not the case 
for most models. 

A practical remedy in most cases is to apply 
the method of successive continuation first intro- 
duced by Doedel, Friedman and Monteiro (1993) 
for point-to-point problems. This method does 
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not guarantee that a connection will be found 
but works well if we start sufficiently close to a 
connection in the parameter space. Here we gen- 
eralize this method to point-to-cycle connections. 

We first consider the case dim IV" = 1. To 
start, we introduce a EVP composed of ([7j), (fT6|) . 
and a modified version of (fT8|) . namely: 



x+ - T+/(x+, a) 


= 0, 


(22a) 


x+(0) -x+(l) 


= 0, 


(22b) 




= 0, 


(22c) 


w + T~^f^{x~^, a)w + Xw 


= 0, 


(22d) 


w{l) - sw{0) 


= 0, 


(22e) 


{w{0),w{0))-l 


= 0, 


(22f) 


u — Tf{u, a) 


= 0, 


(22g) 


(/(x+(0),a),M(l)-x+(0))-/ii 


= 0, 


(22h) 



where \1/ in fl22cl) defines any phase condition fix- 
ing the base point x'^{0) on the cycle O^; for 
example 

=x+(0)-aj, 

where aj is the jth-coordinate of the base point 
at some given parameter values, and hi is a ho- 
motopy parameter. 

Take an initial solution to this EVP that col- 
lects the previously found equilibrium-related da- 
ta, the cycle-related data (x+, T"*") including x"^(0), 
the eigenfunction-related data (w,A), as well as 
the value of hi computed for the initial "connec- 
tion" 

u{t) =e + ewe^"^", r G [0,1], (23) 

which is a solution of the scaled linear approxi- 
mation of ([1]) in the tangent line to the unstable 
manifold W!^ of ^. Ey continuation in (T, hi) for 
a fixed value of a, we try to make hi = 0, while 
m(1) is near the cycle O"*", so that T becomes suf- 
ficiently large. 

After this is accomplished, we introduce an- 



other EVP composed of (^^, ([TBI) , and 





x+ -T+f{x^,a) 


= 0, 


(24a) 




x+(0) -x+(l) 


= 0, 


(24b) 


{wiO),u 


(l)-X+(0))-/l2 


= 0, 


(24c) 


W + T+ 


/J(x+,a)w + Xw 


= 0, 


(24d) 




w{l) - sw{0) 


= 0, 


(24e) 




{w{0),w{0))-l 


= 0, 


(24f) 




u — Tf{u, a) 


= 0, 


(24g) 




a),u{l) - x+(0)) 


= 0, 


(24h) 



where /12 is another homotopy parameter. 

Using the solution obtained in the previous 
step, we can activate one of the system param- 
eters, say ai, and aim to find a solution with 
h2 = hj continuation in (ai,/i2) for fixed T. 
Then we can improve the connection by contin- 
uation in (ai,T), restarting from this latest so- 
lution, in the direction of increasing T. Eventu- 
ally, we fix a sufficiently large value of T and con- 
tinue the (approximate) connecting orbit in two 
systems parameters, say (ai,a2), using the origi- 
nal EVP without any homotopy parameter as de- 
scribed in Section 13. 4[ All these steps are illus- 
trated for the Lorenz example in Section 16. 1[ In 
practice, intermediate continuations in e or other 
system parameters may be necessary to obtain a 
good approximation to the connecting orbit. 

When dimVT" = 2, a minor modification of 
the above homotopy method is required. In this 
case, we replace (fT7|) by the explicit boundary 
conditions 

u{0) - ^ - e{civ^^^ + C2V^^'^) = 0, (25a) 
cl + cl = 1, (25b) 

where e is a small parameter specifying the dis- 
tance between u{0) and ^, v^^^ are two linear- 
independent vectors tangent to of the saddle 
^, and Ci,2 are two new scalar homotopy parame- 
ters. Note that if w = (fi, W2, v^)"^ is a solution to 
iQ with V2 7^ 0, then one can use the normalized 
vectors 













-Vi 








\ I 




\ -V2 1 
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Now consider a BVP composed of ([9]), ( l25i) . and 
(122|) . The initial data for this BVP are the same 
as in the case dimiy" = 1, except for 

Ci = 1, C2 = 0. 

The initial "connection" in this case is 

u{t) = e + ee^^\'^^\ T e [0, 1], (26) 

where A = fu{^,(y), to be used to compute the 
initial value of hi in (122hp . 

By continuation in (T, hi) (and, eventually, in 
(ci,C2,/ii)) for fixed values of all other parame- 
ters, we aim to locate a solution with hi = 0, 
with u{l) near the base point of the cycle so 
that T becomes sufficiently large. We then switch 
to the BVP composed of ([7D, ([25]), and ([21]), and 
we aim to locate a solution with /i2 = 0, by con- 
tinuation in (ci,C2,/i2) for fixed T. When this is 
achieved, we have a solution to the original BVP 
(ED, dm), and (dlD introduced in Section [Ml and 
containing no homotopy parameters. Using this 
BVP, we can continue the approximate connect- 
ing orbit in one system parameter, say ai, with 
T fixed. 

Examples of such successive continuations will 
be given in Section 16.31 where we consider the 
standard model of a 3-level food chain. In that 
section also an alternative BVP formulation for 
( !25|) is given. When one system parameter is var- 
ied, limit points (folds) can be found and then 
continued in two system parameters. 

5 Implementation in AUTO 

Our algorithms have been implemented in auto, 
which solves the boundary value problems using 
superconvergent orthogonal collocation with adap- 
tive meshes, auto can compute paths of solu- 
tions to boundary value problems with integral 
constraints and non-separated boundary condi- 



tions: 

U{r)-F{U{r),(3) = 0, re [0,1], (27a) 
6(f/(0),f/(l),/3) = 0, (27b) 

/ q{U{T),P)dT = 0, (27c) 
Jo 

where 

U{-), F{-, ■) e M"^ b{-, ■) e W"^, q{-, ■) e W% 

and 

Here /3 represents the Ufp free parameters that are 
allowed to vary, where 

nfp = nbc + nic-nd + l. (28) 

The function q can also depend on U and on the 
derivative of U with respect to pseudo-arclength, 
as well as on U, the value of U at the previously 
computed point on the solution family. 

For our primary BVP problem ([7]) or ([9]) with 
f[T6|) or f[T71) . respectively, and f[T8|) . we have 

rid = 9, = 0, 

and = 19 or 18, respectively, since ([7]) and 
are treated as boundary conditions. 

6 Examples 

In this section we illustrate the performance of 
our algorithm by applying it to three model sys- 
tems, namely, the Lorenz equations, an electronic 
circuit model, and a biologically relevant system. 

6.1 The Lorenz system 

One of the best-known dynamical systems that 
has a heteroclinic point-to-cycle connection is the 
three-dimensional Lorenz system, given by 

xi = a{x2-xi), 
±2 = rxi-X2-XiX3, (29) 
^3 = X1X2 - 6x3, 
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Figure 3: Continuation in T: (a) T = 1.43924; (b) T = 1.54543; (c) T = 2.00352. 



with standard parameter values a = 10, b = 
8/3, and where r is the usual bifurcation param- 
eter. With these parameter values, a supercriti- 
cal pitchfork bifurcation from the trivial equilib- 
rium occurs at r = 1, giving rise to two sym- 
metric nontrivial equilibria. At r ^ 13.962 there 
are two symmetry-related orbits of infinite period 
that are homoclinic to the origin, and from which 
two families of saddle cycles arise (together with 
a nontrivial hyperbolic invariant set). A subcriti- 
cal Hopf bifurcation of nontrivial equilibria takes 
place at ^ 24.7368, where these two cycles 
disappear. 

At a critical value Vhet there is a heteroclinic 
point-to-cycle connection, that generates a chaotic 
attractor, see Afraimovich et al. (1977). Its do- 
main of attraction is bounded by the stable invari- 
ant manifolds of the saddle cycles. Beyn (1990) 
found Thet ~ 24.05, and later Died and Rebaza 
(2004) calculated 

rhet = 24.057900322267... 

The heteroclinic connection can be continued 
in two parameters, for example r and a with b 
fixed. The resulting curve in the r, a-plane was 
first shown in Appendix II, written by L.P. Shil'- 
nikov, to the Russian translation of the book by 
Marsden and McCracken (see Pampel (2001), Die- 
ci and Rebaza (2004), for more recent related re- 
sults). As shown by Bykov and Shilnikov (1992), 



the canonical Lorenz attractor appears by cross- 
ing only a part of the heteroclinic connection curve. 

We begin at r = 21.0 and consider a saddle 
limit cycle of ( l29l) with the base point 



X 



(0) = (9.265335,13.196014,15.997250) 



and period T"*" = 0.816222. This cycle can be 
obtained easily by continuation in auto and has 
two nontrivial multipliers: 

/i+ = 0.0000113431, /i+ = 1.26094. 




0.00 0.20 0.40 0.60 

0.10 0.30 0.50 0.70 



Figure 4: Two profiles of the truncated connect- 
ing orbit in the Lorenz system scaled to the unit 
time interval: (a) T = 2.00352; (b) T = 3.0. 

To compute the eigenfunction w, we first con- 
tinue the trivial solution of the BVP fll8al) . fll8b|l . 
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Figure 5: The bifurcation curve of the Lorenz system corresponding to the point-to-cycle connection. 



fITU]) . and fl^Tl) . to detect a branch point at 

A = ln(/i+) = 0.231854, 

from which a nontrivial branch is followed until 
the value /i = 1 is reached. This gives a nontrivial 
eigenfunction with ||w(0)|| = 1, namely, 

w(0) = (0.168148,0.877764,-0.448616)'^. 

In these continuations all problem parameters, 
that is r, a, and 6, are fixed. 

The next step is to find an approximation to 
the connecting orbit. For this, we consider the 
BVP (El), (Iini), and (El with 



=x^(0) - 9.265335 

and continue its solution at fixed system parame- 
ters with respect to (T, h]). Figure [3] shows three 
consecutive solutions with h\ = 0. The end point 
of the last solution (with T = 2.00352) is located 
near the base point x+(0) of the cycle 0+. Us- 
ing this solution as the initial data for the BVP 
([7j), (fT6|) . and (^^, we do a continuation in (r, /i2) 
with T fixed until /12 = is detected. This occurs 
at r = 24.0720, and ensures that the end point of 
the connection is in a plane orthogonal to w{0), 
i.e., in the tangent plane to at x+(0). 



The primary BVP consisting of ([7]), (HM . and 
f[TH]l is used for further continuation runs. First, 
the length of the connecting orbit is increased by 
continuation in (r, T) until T = 3.0. The cor- 
responding parameter value r = 24.0579 gives 
a good approximation for r^et, since the 'tail' of 
the connecting orbit follows the cycle several 
times; (see Figure H]). 

Finally, continuation in the two system pa- 
rameters (r, a) with T fixed, gives the bifurcation 
curve corresponding to the point-to-cycle connec- 
tion in fl2^ . see Figure O 

6.2 A circuit model 

The next example is one from the HOMCONT de- 
mos of Champneys et al. (1999), namely, the elec- 
tronic circuit model of Freire et al. (1993; see also 
the AUTO demos tor and cir). The equations are 



rxi = -{/3+iy)xi+/3x2- 
±2 = pxi - {P+'j)x2 - X3 
is = X2, 

(30) 

where 7 = 0, r = 0.6, 03 = 0.328578, 63 = 0.933578, 
and u and (3 are bifurcation parameters. With 
HoMCONT it was shown previously that a homo- 
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clinic connection to the origin occurs for 



1.5_ 



Uinit = -0.721309 , (3init = 0.6 

with truncated time interval T = 200. Continua- 
tion in two-parameter dimension then leads to a 
Shil'nikov-Hopf bifurcation at 

u = -1.026445 , p = -2.330391 ■ 10'^, 

where a limit cycle bifurcates from the equilib- 
rium, effectively turning the homoclinic connec- 
tion into a heteroclinic one (see auto demo cir). 
We can now compare the results from the con- 
tinuation in HOMCONT with the results from the 
application of our BVP system. 

The equilibrium in this system is a saddle- 
focus, and we therefore have n~ = 2 and n~ = 1. 
To generate appropriate starting data we locate a 
Hopf bifurcation, with (3 as free parameter, from 
where a cycle is continued up to a selected value 
of /3, say, (3 = —0.32. The saddle limit cycle O"*" 
has the base point 

x+(0) = (0.03448278,0.46460323,0.4737975) 

and period T"*" = 6.3646138. The nontrivial mul- 
tipliers are 

/i+ = 3.986051 ■ 10^^ /i+ = 18.85438 

The eigenfunction of this cycle is computed as 
described in Section 14. 2[ which yields 

w{0) = (0.99950,-0.019205,0.024767)^ 

and the log multiplier 

A = -13.579343187. 

An approximation of the connecting orbit is 
then obtained using BVP dB, (16]), and ([22]), with 

^[a;+] =x+(0) - 0.46460323. 

The software content is used to get a good 
approximation of the connection period T, after 
which shooting in MATLAB is used to obtain the 
orbit itself for the given period. 



X2 




-3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 

Xi 

Figure 6: A point-to-cycle connection of the elec- 
tronic circuit model, projected onto the Xi,X2- 
plane. 

Continuation of this approximate orbit with 
respect to {T,hi) yields several orbits with hi = 
0. For T = 11.59816 the orbit is close enough 
to the X2 base coordinate to use the data for the 
BVP dZ]), (dS]), and Continuation in (i^, /is) 

is done until a zero of h2 is reached. 

The primary BVP ([7]), ([I6]), and ([H]) is used 
in the subsequent computations. Continuation in 
{u, T) gives orbits of any desired period T; we 
used T = 20 with 

u = -1.500498. 

At this point continuation can be done in (z/, /3). 

In Figure [6] we see a point-to-cycle connection 
in a xi,a;2-plot at some selected parameter values. 
It is apparent that the homotopy method has re- 
sulted in a good approximation of the connect- 
ing orbit. Figure [7] shows the composite results 
of the two-parameter continuation of the homo- 
clinic connection in HOMCONT and our continu- 
ation of the heteroclinic connection. Label 5 is 
the starting point of the continuation of the ho- 
moclinic connection that terminates at the solu- 
tion labelled 1. Beyond this solution Homcont 
gives spurious results. Note that label 1 coincides 
with label 9, where the curve of the heteroclinic 
connection turns back onto itself, i.e., the con- 
tinuation reverses direction approximately at the 
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Figure 7: Continuation in (z/, j3) of the point- 
to-cycle connection, as explained in detail in the 
text. 

point where the Shil'nikov-Hopf bifurcation oc- 
curs. Plots in AUTO of the limit cycle data (not 
shown) reveal that indeed the cycle shrinks prac- 
tically to a point, before the continuation reverses 
direction. 

6.3 A food chain model 

The following three-level food chain model from 
theoretical biology is based on the Rosenzweig- 
MacArthur (1963) prey-predator model. The equa- 
tions are given by 

Xi = Xi(l - Xi) - /i(xi,X2), 
2^2 = fl{Xl,X2) - f2{x2,X3) - diX2, (31) 
is = f2{x2,X3) - d2X3, 

with Rolling Type-II functional responses 



1 + biU 



1,2. 



The death rates di and d2 are used as bifur- 
cation parameters, with the other parameters set 
at ai = 5, 02 = 0.1, 6i = 3, and 62 = 2. 

It is well known that this model displays chao- 
tic behaviour in a given parameter range, see Ho- 
geweg and Hesper (1978), Klebanoff and Hastings 
(1994), McCann and Yodzis (1995), Kuznetsov 
and Rinaldi (1996), and Kuznetsov et al. (2001). 



Previous work by Boer et al. (1999, 2001) 
has also shown that the regions of chaos are in- 
tersected by homoclinic and heteroclinic global 
connections. In particular, a heteroclinic point- 
to-cycle orbit connecting a saddle with a two- 
dimensional unstable manifold to a saddle cycle 
with a two-dimensional stable manifold can exist. 
It was shown that the stable manifold of this limit 
cycle forms the basin boundary of the interior at- 
tractor and that the boundary has a complicated 
structure, especially near the equilibrium, when 
the heteroclinic orbit is present. These and other 
results were obtained numerically using multiple 
shooting. In this section we reproduce these re- 
sults for the heteroclinic point-to-cycle connec- 
tion. Using our homotopy method we obtain an 
accurate approximation of the heteroclinic orbit. 
A one-parameter bifurcation diagram then shows 
limit points, which correspond to tangencies of 
the above-mentioned two-dimensional manifolds. 
We then continue the limit points in two param- 
eters. 

A starting point can be found, for example, at 
di ^ 0.2080452, 4 = 0.0125, where there is a fold 
bifurcation in which two limit cycles appear. This 
also corresponds to the birth of the heteroclinic 
point-to-cycle connection. 

Before using the homotopy method to obtain 
an approximation of the point-to-cycle connec- 
tion, we locate a Hopf bifurcation, for instance 
at di ^ 0.51227, ^2 = 0.0125. The limit cycle 
born at this Hopf bifurcation is continued up to 
a selected value of di, say, di = 0.25. 

We now have an equilibrium 

C = (0.74158162, 0.16666666, 11.997732) 

and a saddle limit cycle with the base point 

x+(0) = (0.839705,0.125349,10.55289) 

and period = 24.282248. Its nontrivial multi- 
pliers are 

/i+ = 0.6440615, /i+ = 6.107464 ■ lOl 

The eigenfunction w is obtained as described 
in the previous sections. Continuation of the triv- 
ial solution of the BVP ffTSil) . ffTSb]) . (fT9|). and 
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( 12T|) and the subsequent continuation of the bi- 
furcating family until h = 1, yields the multiplier 



A = ln(/i^ 



-0.439961. 



Note that we use the stable multiplier, because of 
the projection boundary conditions. The associ- 
ated nontrivial eigenfunction w{t) with ||tf;(0)|| = 
1 has 

w{0) = (0.09306,-0.87791,-4.69689)^. 

We now consider a B VP composed of ([9]) , (l25h) , 
and (l22l) . Using content and matlab we ob- 
tain an approximation of the connection with the 
boundary condition 

'^[x'^] = x^(0) — 0.125349 Figure 8: An approximation to the point-to-cycle 

connection projected onto the (x2, X3)-plane for 
and period T = 155.905. The starting point is ^^e food chain model with = 5, = 0.1, h = 
calculated by splitting the normalized adjoint sta- 3 ^2 = 2 (ii = 25 and d2 = 0125 
ble vector (evaluated at di = 0.25, ^2 = 0.0125) 

V = (0.098440,0.168771,0.0049532)^ 




3:3 



into v^^^ and v^'^\ as described in Section H73| and 
multiplying it by a small e, say e = 0.001. In our 
case the starting point was 

m(0) = (0.742445,0.166163,11.997732). 

The first homotopy step involves continuation 
in {hi,T). However, this does not lead to zeroes 
of hi. To obtain hi = we expand the previous 
set of B VPs with (p5b). Subsequent continuation 
in (ci,C2,/ii) gives a solution with hi = that 
indeed ends near the base point x^(0) of the limit 
cycle. 

For continuation in the second homotopy step, 
a switch is made to a BVP composed of Qj, ([21]) 
and fl2S]) . Continuation in (ci,C2,/i2) leads to 
some solutions with /i2 = 0. 

The obtained approximate connecting point- Figure 9: Several point-to-cycle connections in 
to-cycle connection now suffices for continuation the food chain model with different values of di. 
in system parameters. Before doing a continua- 
tion in a system parameter the connection is im- 
proved by increasing the connection period. A 
user-defined point of T = 300 suffices. Next, 




0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 

2^2 
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the parameter e is decreased up to a user-defined 
point oi e = —1 ■ 10~^, so that the starting point 
m(0) is shghtly away from the equilibrium ^. Fig- 
ure [H] displays a projection of the point-to-cycle 
connection onto the (x2, a;3)-plane. 

Now the connecting orbit can be continued up 
to a limit point in one system parameter. Fig- 
ure [9] displays three connecting orbits obtained af- 
ter continuation with respect to ai = di. Contin- 
uations in di result in the detection of the points 

di = 0.280913 and di = 0.208045 

where the first one is a limit point and the sec- 
ond one a termination point. This point coincides 
with a tangent bifurcation for the limit cycle to 
which the point-to-cycle orbit connects. Contin- 
uations in d2 result in the detection of the points 

d2 = 0.0130272 and rfa = 9.51660 ■ 10"^ 

which are both limit points. Any of the detected 
limit points can now be used as a starting point 
for a two-parameter continuation in a = di or d2- 
In practice, the connection period may have to 
be increased or decreased to obtain the full two- 
parameter continuation curve. In the demo, the 
last limit point {d2 = 9.51660-10"^) is the one se- 
lected for the food chain model. The two-parame- 
ter continuation curve terminates at both ends 
in codim 2 points lying on the above-mentioned 
tangent bifurcation for the limit cycle. These 
points coincide with the log multiplier A = 0. 
Observe that this corresponds to the point di = 
0.208045, detected in the one-parameter continu- 
ation, where also A = 0. 

For the continuation in two system parame- 
ters, the BC ( 12^ proves ineffective, since it leads 
to the detection of several spurious limit points. 
This is, because the orbit spiralling out from the 
equilibrium has an elliptical shape. The circle of 
a small radius, centered at the equilibrium, inter- 
sects the spiral at several points, one of which is 
the starting point of the connecting orbit. During 
continuation, with a changing problem parame- 
ter, the spiral will change size and the starting 



point on the circle may collide with another such 
point where the circle and the spiral intersect. 
This intersection would correspond to a fold with 
respect to the problem parameter. As a result, to 
obtain a full continuation curve of the connecting 
orbit in two system parameters, some restarts are 
required. 

In order to avoid these spurious folds, we re- 
turned to the original BC (1171) with (117bp in the 
form 



(32) 



where j is either 1,2 or 3. By setting j = 2 we are 
in line with the work by Boer et al. (2001), who 
used a Poincare plane through the equilibrium ^ 
where X2 = C,2- 

Figure [10] shows the curve of limit points T/jgt 
that is computed with the method described above, 
using the standard switching and fold-following 
facilities of AUTO. This curve can be obtained 
in one run, given the connection period is chosen 
conveniently. It agrees with the results previously 
obtained by Boer et al. (1999) by labourious mul- 
tiple shooting. 

7 Discussion 

Our continuation method for point-to-cycle con- 
nections, using homotopies in a boundary value 
setting is both robust and time-efficient. Detailed 
AUTO demos that carry out the computations de- 
scribed in Section 6 are freely downloadable from 
www . bio . vu . nl/thb/ research/ pro j ect/ globif . 

Although the method was presented for 3D- 
systems, it can be extended directly to point-to- 
cycle connections in n-dimensional systems, when 
the unstable invariant manifold of the equilibrium 
^ is either one-dimensional or has codimension 
one, while the stable invariant manifold of the 
cycle 0+ has codimension one. 

In the forthcoming Part II of this paper, we 
will extend our method to include detection and 
continuation of cycle-to-cycle connections. 
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Figure 10: A two-parameter bifurcation diagram of the food chain model that shows the region where 
there exist point-to-cycle connections. The region is bounded on one side by the cycle fold, (Tc), and 
on the other side by the curve Thet, the locus of limit points of the heteroclinic connections. 
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A ]V[onodrorQy niclt rices Together with flMl) . consider the adjoint vari- 
ational system 

In order to approximate the invariant manifolds . ^ ^ 

of a hmit cycle we use eigenvalues and eigenfunc- ^ ~ ~^ ^^^^ ' C ^ ^ (37) 

tions of appropriate variational problems. These , ^ • • -i- i i u 

, ^ ^ ^ and the correspondmg matrix mitial-value prob- 
eigenvalues m turn are the eigenvalues of the so- 



Z = -A^{t)Z , ZiO) = /„, (3^ 



called monodromy matrix. 

To define an eigenfunction of the periodic so- 
lution x{t+T) = x{t), where T is the period of the which is the adjoint system to Note, that 
cycle, of an autonomous system of smooth ODE's the multipliers of the adjoint monodromy matrix 

u = f{u), (33) N = Z{T) 

write a solution of this system near the cycle in ^^^^^^^ multipliers of the monodromy ma- 

the form ^^^^ ^ ~ ^(^)- The proof of this well-known fact 



u{t) = x{t)+m 



goes as follows. Compute 



where ^(t) is a small deviation from the periodic d /^Ty-v _ y -\- Z^^^ 

solution. After substitution and truncation of dt dt dt 

the 0(||^|p)-terms, we obtain the following vari- = [—A^Z)^Y + Z^ AY 

ational system: _ z'^{—A)Y + Z^AY = 

^ = A{t)^ , e G , (34) Since Z{0) = Y{0) = 4, we get Z'^{T)Y{T) = 4, 

which implies 

where A(t) = fu{x{t)) is the Jacobian matrix 

evaluated along the periodic solution; A{t + T) = [M"^]""". 

= ^ 

Now, consider the matrix initial- value prob- Due to (IMD, a multiplier /i satisfies v{T) = 

\qjj^ Att'(O) with v{0) or, equivalently, it is a solu- 

Y = A(t)Y , ^(0) = In, (35) tion component of the following BVP on the unit 

interval [0, 1]: 

where /„, is the unit n x n matrix. Its solution 

Y{t) at t = T is the monodromy matrix of the ( it -\- TA{t)v = , 

cycle: I v{l) - fiviO) = , (39) 

M = Y{T). \ (t;(0),t;(0))-l = 0. 

The monodromy matrix is nonsingular. Any so- p^^^^ ^^^^^^ ^^lat /i > and write 
lution ^{t) to flM|) satisfies 



^(T) = M^(0) . (36) 



/i = e^, v{t) =e^^w{t). 
Then w satisfies a periodic BVP, namely: 



The eigenvalues of the monodromy matrix M are 

called the Floquet multipliers of the cycle. There j + TA{t)w + \w = , 

is always a multiplier +1. Moreover, the product 

{ w{l)-w{0) = 0, (40) 

of all multipliers is positive: ( {w{0),w{0))-l = 0. 

r-r \ Similarly, when /i < 0, we can introduce 



/ii/i2 ■ ■ ■ f^n = exp div f{x(t)) dt^ 



^ = -e^, v{t) = e^^w{t) 
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and obtain an anti-periodic BVP 

w + TA{t)w + Xw = , 

w{l) + w{0) = 0, (41) 
{w{0),w{0)) - 1 = 0. 

This technique can easily be adapted to the mul- 
tiphers of the adjoint variational problem (!37|) . 

Finally, we note that the eigenvalue problem 
for a Floquet multiplier 

Mv - fiv = 

can be considered continuation problem with 
n + 1 variables {v, fi) G M" x M defined by n equa- 
tions. This continuation problem has a trivial 
solution family = (0,yu). An eigenvalue 

/ii corresponds to a branch point, from which a 
secondary solution family {v,fii) with v ^ em- 
anates. This nontrivial family can be continued 
using an extended continuation problem 

( Mv - fiv = , 
\ {v,v)-h = 0, 

which consists of n + 1 equation with n + 2 vari- 
ables h). If /i = 1 is reached, we get a nor- 
malized eigenvector v corresponding to the eigen- 
value /ii, since along this branch = /ii. Gen- 
eralization of this procedure to the BVP fl39|) (as 
well as to fHOj) . pTl) . and their adjoint versions) 
is straightforward. 
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